The sensitivity of the vortex filament method to different reconnection models 
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We present a detailed analysis on the effect of using different algorithms to model the reconnection 
of vortices in quantum turbulence, using the thin-filament approach. We examine differences between 
four main algorithms for the case of turbulence driven by a counterflow. In calculating the velocity 
field we use both the local induction approximation (LIA) and the full Biot-Savart integral. We 
show that results of Biot-Savart simulations are not sensitive to the particular reconnection method 
used, but LIA results are. 

(N ; 

^ : 

C/5 ! I. INTRODUCTION 

i> : 

Turbulence in the quantum, low temperature phase of liquid helium ( 4 He), also known as quantum turbulence 
[ll-Q, consists of reconnecting quantized vortex filaments, arranged in a random, disordered tangle. Due to quantum 
mechanical constraints, each vortex carries the same fixed circulation h/m, where h is Plack's constant and m is the 
mass of one atom. This quantity is called the quantum of circulation k. A number of experimental methods have been 
used to create this form of turbulence, for example agitating superfluid liquid helium with propellers [H, Hj], forks @, 
M ■ or grids |7j; these techniques are also used to create turbulence in ordinary fluids. 

At finite temperatures superfluid helium is a two fluid system: a viscous normal fluid component coexisting with 
an inviscid superfluid component. The superfluid vortices interact with the thermal excitations which make up the 
' normal fluid, thus introducing a mutual friction force between the two fluid components. This means that turbulence 
. in the quantum fluid can be driven by the flow of the normal fluid or vice-versa. One particular example is superfluid 
turbulence driven by a heat flow [sMlOjj . Note that this form of turbulence has no classical anal ogy . 

Recently new flow visualization techniques, such as tracer particles [ll|,[l|j], Andreev scattering [is| and laser-induced 
I ■ fluorescence [l4j, have added substantially to our knowledge of the nature of quantum turbulence. Even with these 
advances in flow visualization numerical simulations will continue to play a crucial role in furthering our knowledge. 
Indeed, building on the pioneering work of Schwarz [l5j . numerical simulations have always been important in the 
fiel d ll6H26l | , and the recent experimental progress has highlighted their importance in interpreting experimental data 
(27l |28| |. This article is concerned with verifying that the vortex reconnection procedure used in the popular vortex 
filament method (VFM) is robust. Before we discuss the modelling of reconnections in the VFM, we shall give a brief 
outline of the VFM method. 

' In superfluid helium, the vortex core radius («o ~ 10~ 8 cm) is many orders of magnitude smaller than the average 
separation between vortex lines (typically from 10~ 2 to 10~ 4 cm) or any other relevant length scale in the flow. Starting 
from this key observation, Schwarz [15| modelled vortex lines as spaces curves s = s(£ i i) of infinitesimal thickness, 
where t is the time and £ is arc length, using the classical theory of vortex filaments [29[. In the VFM these space 
. curves are numerically discretized by a large, variable number of points Sj (i — 1, • • • N), which hereafter we refer to 
as vortex points. In all implementations of the VFM to supcrfluids the number of vortex points varies in time, as it is 
desirable to maintain a relatively constant resolution along the filaments by adding or removing points as the length 
of the filaments changes. In a recent paper [25[ we have discussed the various numerical details at length, here it is 
J> \ suffice to say that the resolution along the filaments lies between an upper, 6, and lower, 5/2, bound. 

The motion of a vortex filament is determined by the normal fluid, through mutual friction, and by the induced 
velocity of the other vortices. The governing equation of motion of the superfluid vortex lines, at point s is given by 
■ the Schwarz equation [lj| 
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— = v s + as' x (v w - v s ) - a's' x [s' x (v n - v s )] , (1) 
at 

where a, a' are temperature dependent friction coefficients [13,[3l|, v„ is the normal fluid's velocity, and the velocity 
v s is the velocity field the quantised vortices alone induce. The prime denotes derivative with respect to arclength, 
e.g. s' = ds/dt;. Note that, strictly speaking, one should also model the back reaction of the superfluid on the normal 
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FIG. 1: Schematic reconnection procedure for Type I and II reconnections described in section (TT] If two vortex points are 
closer than a critical distance A = 8/2 then a flag swapping operation changes the topology of the filaments and performs the 
reconnection. From (1 — > 2 — > 3) and (4 — > 5 — > 6) to (6 — > 5 — > 2 -> 3) and (1 -> 4). 



component; v n should be the solution of the Navier-Stokes equation modified by mutual friction. However, for the 
sake of simplicity, here we shall simply prescribe the normal flow. 

In our system the quantised vortices define the vorticity field; we recover the velocity field v s by numerically solving 
the Biot-Savart (BS) integral, 



4tt J c |s - 



x dr. 



(2) 



where the line integral extends over the entire vortex configuration C. We de-singularize the integral in a standard 
way (l5j . 

Schwarz [lj| proposed the use of the Local Induction Approximation (LIA) [32], [33| as an alternative to the Biot- 
Savart method, which is computationally very expensive. The LIA ignores the non-local contribution to the motion 
of a section of the vortex filament. Instead a vortex line at the point Si moves along the binormal to Sj, with the 
following equation of motion 



dt 
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where c is a constant of order unity and R is the local radius of curvature |s"| _1 of the vortex filament. We shall test 
various reconnection algorithms using both the full BS law, and the LIA in a series of numerical experiments. Again, 
for further details about spatial derivatives and methods to time-ste p th e vortices, we refer the reader to pBj . 

We know from experiments [34j and from more microscopic models j35l - [38l | that superfluid vortex lines can reconnect 
with each other when they come sufficiently close, as envisaged by Feynman [39j . Superfluid vortex reconnections 
do not violate Kelvin's theorem as near the axis of the vortex core, where density and pressure vanish and velocity 
diverges, the governing Gross-Pitaevski equation (GPE) differs from the classical Euler equation. 

Whilst vortex reconnections are natural solutions of the GPE, within the VFM reconnections must be modelle d by 
supplementing Eq. @ with an algorithmical reconnection procedure. This was originally proposed by Schwarz [l5| . 
and since then a number of alternative algorithms have been proposed. Until now no detailed test of the effects of 
varying this procedure has been performed. Whilst a number of studies have shown good agreement between results 
using the VFM and experimental results [26|, |4fJ, |4l| , a detailed study of the reconnection algorithm is timely; that is 
the purpose of this study. 



II. RECONNECTION ALGORITHMS 

In his pioneering paper [i~o| . Schwarz's suggested that vortex lines reconnect whenever the distance between a pair 
of vortices is less than A = 2R/[cln(R/ao)], where R is the radius of curvature at the reconnection point, and c is a 
constant of order unity. However, this approach can lead to non-physical reconnections. Consider, for example, two 
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FIG. 2: Schematic reconnection procedure for Type III reconnection algorithm described in section[II] The segments (1 — ► 2 - 
3) and (4 — > 5 — ► 6) evolve into (6 — ¥ 3) and (1 — > 4). Points 2 and 5 are eliminated. 





FIG. 3: Schematic representation of the motion of vortex filaments under the Type IV algorithm. If the segments will collide 
(as shown) then the reconnection is performed as in Fig. [T] 



almost straight vortices; under Schwarz's criterion these vortices must reconnect even if they are very far apart, as a 
large radius of curvature R, results in a large value for A. In this work we therefore avoid the use of this reconnection 
criterion, instead focusing on methods which have been used in more recent studies. 

Most recent studies have related the reconnection distance to the space resolution along the filaments [HI, |43| . We 
consider three differing reconnection algorithms which take the critical reconnection distance to be related to the 
maximum spatial resolution, A = 5/2. We define the first reconnection algorithm as Type I. In a Type I reconnection 
vortices are simply reconnected if their separation if less than A. Motivated by the fact that a reconnection is 
dissipative event, leading to phonon emission 4J|, we consider an improved algorithm denoted Type II. As vortex line 
length is a proxy for the kinetic energy, for a Type II reconnection not only must the distance between the reconnecting 
filament be less that 5/2, but also the line length must be reduced by the change of topology. A schematic for the 
Type I and II algorithms can be seen in Fig. Q] The third reconnection model, Type III, can be considered an 'ultra' 
dissipative algorithm. Type III differs from Type I because the points which triggered the reconnection are eliminated 
leading to a greater loss of line length, and hence energy. A schematic of this algorithm is displayed in Fig. [5] 

In a recent work Kondaurova and Nemirovskii 21| introduced a new model (Type IV) of reconnections within 
the VFM, which tests whether line segments will meet during the time-step. Candidate pairs of vortex points (i & 
j), points which are close enough to reconnect, are identified. The Type IV algorithm then assumes that the line 
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segments between each of the pair will move with a constant velocity (v(s,) & v(sj)) during the time-step. Under 
this assumption a set of simultaneous equations can be constructed, 

Xi + v x (si)At + (x i+ i - Xi)ip = Xj + v x (sj)At + (xj+i - Xj)4> 
Vi + v y (si)At + (y i+1 - yi)i> = y. } + v y (sj)At + (y j+1 - y^cj) 
Zi + v z (s l )At + (zi+i - Zi)ip = zj + v z (sj)At + (zj+i - Zj)4> 



O<0<1,O<-0<1 



such that if a solution for cj> and ip can be found, then the filaments will collide during the time-step. Here s, = 
(Xi,yi,Zi), s i+ i = (xi+i,yi+i,Zi+i) Sj = (xj,yj,Zj) and Sj+i = (x j+ i,yj + i,Zj + i) are the coordinates of the pair of 
points and their neighbours along the filament, and At is the numerical time-step. If the line segments will meet, then 
the filaments are reconnected in the same manner as in Fig. [TJ We solve for i/> and <j> algebraically using the equations 
for the x and y components. The values of tp and <fi obtained are then tested in the equation for the z component by 
verifying that, 

\zi + v z (si)At + (zi+i - Zi)tp - zj + v z (sj)At + (zj+i - Zj)4>\ < e, (4) 

where e is the imposed the tolerance of the method. We define candidate points (i & j) as vortex points within 
distance A = S. A schematic diagram of the Type IV algorithm is shown in Fig. [3] 

Before performing any reconnection we test the distance from a point i to all other points, which are not the nearest 
neighbours along the filament. We then begin the test for a reconnection based on the closest reconnection point, 
assuming the distance is smaller than A. This means self-reconnections (which can arise if a vortex filament has 
twisted by a large amount) are treated in the same manner as reconnections between different filaments. Finally as 
reconnections must preserve the orientation of the vorticity, we check that the two reconnecting filaments are not 
parallel. To do this we form local (unit) tangent vectors, for example if the reconnecting points are i and j we 
can readily calculate and s^-. We then test that • < 0.965, which ensures that the minimum angle between 
reconnecting filaments is approximately 15 degrees. We have tested that altering this threshold to smaller angles does 
not lead to any discernible difference in the results, hence for brevity all results here are subject to this criterion. 

III. NUMERICAL SIMULATIONS 

We now the numerical simulations used to test the various reconnection algorithms described in the previous 
section. As a benchmark we choose counterflow turbulence, the relative motion of the normal fluid and superfluid 
components sustained by an applied heat flow in the direction of v„. The superfluid flows in the opposite direction 
so that p n v n + p s v s = 0, where p n and p s are respectively the normal fluid and superfluid densities. We choose 
this form of turbulence as there has been a wealth of experimental studies [1, [13, [HI , as well as some recent detailed 
numerical simulations [40l . [4(| . Counterflow turbulence was also used in a recent numerical study by Kondaurova and 
Nemirovskii [13| , where they reported that the use of a reconnection method, similar to the Type IV algorithm, gave 
a steady state solution for LIA simulations. We shall also test these claims in this study. 

Our calculations are performed in a periodic cube with sides of length D = 0.1 cm. Superfluid and normal 
fluid velocities v„ and v s are imposed in the positive and negative x directions respectively, where v ns = |v„ — v s | is 
proportional to the applied heat flux. Simulations are performed with the same numerical resolution S = 1.6 x 10~ 3 cm 
and time-step At — 10~ 4 s. The initial condition consists of a set of vortex loops set at random locations, which are 
displayed in Fig. [4j 

During the evolution we monitor the vortex line density calculated as, 

L=±, (5) 

where V = D 3 is the volume of the computational domain, C is the entire vortex configuration and A = * s ^ ne 

total vortex line length. We also monitor the number of reconnections per unit time (the reconnection rate), £. 

We run simulations, with the different reconnection algorithms, using both the LIA and the full BS law with 
v ns = 0.55 cm/s and temperature T = 1.6 K (a = 0.098, a' — 0.016). We then run further simulations, using only the 
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FIG. 5: The vortex line density L (cm) plotted as a function of time t (cm) for simulations with the LIA (left) and the full 
BS law (right) with different reconnection algorithms, v ns = 0.55cm/s and T = 1.6K. The symbols are as follows: (magenta) 
diamonds: Type I reconnection; (yellow) triangles: Type II; (cyan) pentagons: Type III; (black) circles: Type IV (e = 10 -4 ); 
(blue) squares: Type IV (e = 10 -3 ); (red) crosses: Type IV (e = 10~ 2 ) and (green) asterisks: Type IV (e = 10 _1 ). 



full BS law, for three further values of v ns . This allows us to test the theoretically predicted Q and experimentally 
and numerically (40l | verified law for the steady state line length, 

L = l 2 v 2 ns , (6) 

where 7 is a temperature-dependent parameter. Our aim is to determine whether the different reconnection algorithms 
yield different values of 7. We stress that 7 is a macroscopic quantity which is known from experiments. 
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FIG. 6: The reconnection rate, £ (s _1 ), scaled by the vortex line length, A (cm), plotted as a function of time, t (s), for the 
LIA (left) and BS (right) simulations; plotting symbols and colors is as in Fig. [S] Note that in the case of the LIA simulations 
a number of the simulations reach a degenerate steady state, as is visible in Fig. [5] (left). In this case the reconnection rate 
drops to zero and hence the results are not visible on this semi-log scale plot. This does not happen in the full BS simulations. 



IV. RESULTS 



The results for the initial simulations using both the LIA and the full BS law are presented in Figs. [SH5] Figure 
[5] shows the vortex line density plotted as a function of time for the LIA (left) and the BS (right) simulation, using 
various reconnection algorithms. The results for the Biot-Savart law are very encouraging showing little difference 
between the simulations with the different reconnection algorithms. In all simulations there is an initial rapid growth 
in the vortex line density, followed by a more gradual increase until eventual saturation to a fluctuating steady state. 

As we discussed earlier, the reconnection method is the one ad-hoc aspect of the VFM, and it seems that, at least 
for quantum turbulence driven by counterflow, the results are very robust. Note infact that all calculations performed 
with the BS law (Fig.[5]right) converge to approximately the same value of L. We believe this is a very useful result as 
this means that one can confidently draw conclusions from simulations which use the BS method. We note, however, 
that one should probably always err on the side of caution with VFM simulations and check the effect of different 
reconnection algorithms with a different numerical experiment. 

Figure [6] (right) shows the number of reconnections per unit time per unit length, £/A, for the BS simulations, and 
the picture is the same. For each of the different algorithms the reconnection rate per unit line length is approximately 
the same. The only exception is the simulation with the Type I reconnection algorithm, which shows a slightly higher 
reconnection rate. This is understandable as this is the algorithm which has the least restrictions and simply reconnects 
filaments if the distance between them is less than 5/2; however there is no discernible impact on the evolution of the 
vortex line density. Snapshots of the system at the end of simulations are plotted in Fig. [3] We shall return to discuss 
the further simulations using the BS method at the end of this section. 

The picture for the LIA simulations is very different. We find that even with the same initial conditions, slight 
differences in the reconnection algorithms lead to large changes in the overall evolution of the system. Kondaurova 
and Nemirovskii [47j | found that using a reconnection scheme similar to the Type IV algorithm used here they reached 
a steady state. We find that certain simulations do reach a steady state, although the steady states we find are 
degenerate. By degenerate we mean that the system evolves to a set of very straight vortices, arranged in planes 
parallel to the counterflow direction. These planes are then simply advected in the direction of the counterflow with a 
negligible change in vortex line density. In this configuration the reconnection rate drops to zero. This is the eventual 
fate of the simulations with Type I, II and IV (e = I CP 3 , 10~ 2 ) algorithms. 

In the other simulations a steady state is not found and large fluctuations in the line density are seen. We find the 
same bundles of vortex lines, stratified in layers, observed by Schwarz [l5[ and recently reproduced by Adachi [io| . 
Schwarz recognised that this is a spurious effect and introduced his artificial mixing procedure to avoid it [To| . see 
Fig. [7] (right). The reconnection rate per unit line length, Fig. |H] (left), also show that LIA simulations are highly 
sensitive to the reconnection algorithm used. We thus agree with other authors that, due to the absence of the 
non-local component of the velocity field, LIA is unsuitable for simulations of turbulence [4C|. 

We now return to the BS simulations, focusing on the affect of varying the counterflow velocity v ns . We are 
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FIG. 7: Snapshots of the vortex filaments from the simulations using the LIA, as plotted in Fig. [5] (left). The left panel shows 
the typical state of each of the simulations at the peak of vortex line density, at t ~ 6(s); this particular plot is taken from the 
simulation using the Type III reconnection algorithm. The middle panel shows the degenerate solution taken at t — 80 s for the 
simulation with Type IV reconnection algorithm (e = 10 -3 ). The right panel shows a plot in the plane perpendicular to the 
counterflow for the simulation with the Type IV (e = 10 _1 ) algorithm at t — 60s: this layered vortex structure was recognised 
by Schwarz [l5l ] and Adachi et al. [iol ]. 




FIG. 8: Snapshots of the vortex filaments in the plane perpendicular to the counterflow for BS simulations, from Fig. [5] (right), 
each at t = 10 s. The left panel is for the simulation with the Type II reconnection algorithm, Type III (middle) and Type IV 
(e = I)" 4 (right). 



interested in the relationship between the steady state line density and the counterflow velocity. As above, the 
temperature is fixed at 1.6K Figure [5] (right) is very representative of the time series plot for the line density, at 
different temperatures, and so we do not reproduce further plots for each of the different values of v ns . Instead in 
Fig. [9] we display the linear relationship between y/~L and v ns , according to Eq. [5] Again, we find that the result is 
not sensitive to which reconnection algorithm is used. 

Finally it is instructive to compare the effect of the reconnection algorithm used on the scaling parameter, 7, and 
compare with experimental results. Our aim here is not to claim that the algorithm which yields the value closest 
to that experimental value should be used in the VFM. Instead we wish to add further weight to our claims that 
results from the VFM can be informative, so long as a physically justifiable reconnection algorithm is used. Table 1. 
shows the values of 7 for the different reconnection algorithms. The experimental value of 7 at this temperature is 
93 0, Ell ; Adachi et al. [Io| gave numerical results of 7 = 109.6. 
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FIG. 9: A plot of the root of the vortex line density, L, and the magnitude of the counterflow velocity, v ns . The linear 
relationship between \[L and v„ a for each of the reconnection schemes is apparent, plotting symbols and colors is as in Fig. [5] 
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Type I 




100.6 




Type II 




111.5 




Type III 


111.7 


Type 


IV (e = 


10" 4 ) 


104.1 


Type 


IV (e = 


lO" 3 ) 


106.8 


Type 


IV (e = 


lO" 2 ) 


111.6 


Type 


IV (e = 


10" 1 ) 


116.4 



TABLE I: A table of values of 7 for the different reconnection algorithms in simulations using the Biot-Savart law at T=1.6K 



V. CONCLUSIONS 

In conclusion we have shown that the vortex filament method is very robust to the reconnection algorithm used, 
provided one uses the full Biot-Savart integral to determine the velocity field. In contrast results from simulations 
using the local induction approximation are highly dependent on the algorithm used, lending further weight to the 
criticism of this method. 
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